function [col, row] = project(rpc_model, x, y, z)
%hsRPC.project 应用RPC模型进行投影, 将物方/像方坐标转换为像方/物方坐标
%
% Syntax:
%   [col, row] = hsRPC.project(rpc_model, x, y, z);
%
% Params:
%   - rpc_model [required] [struct] RPC结构体, 由 hsRPC.read 或者 hsRPC.fitting 获取
%   - x         [required] [numeric] 物方坐标的经度 / 像方坐标的列
%   - y         [required] [numeric] 物方坐标的纬度 / 像方坐标的行
%   - z         [required] [numeric] 物方坐标的高度
%
% Return:
%   - col 像方坐标的列(0-起点) / 物方坐标的经度
%   - row 像方坐标的行(0-起点) / 物方坐标的高度
%
% Author: 
%   iam002, 2025年2月28日
%
% Update: 
%   iam002, 2025年2月28日
%
% Matlab Version: R2023a

    x = (x - rpc_model.longOffset) / rpc_model.longScale;
    y = (y - rpc_model.latOffset) / rpc_model.latScale;
    z = (z - rpc_model.heightOffset) / rpc_model.heightScale;

    num = auxCompute(rpc_model.lineNumCoef, x, y, z);
    den = auxCompute(rpc_model.lineDenCoef, x, y, z);
    row = num ./ den;
    row = row.*rpc_model.lineScale + rpc_model.lineOffset;
    
    num = auxCompute(rpc_model.sampNumCoef, x, y, z);
    den = auxCompute(rpc_model.sampDenCoef, x, y, z);
    col = num ./ den;
    col = col.*rpc_model.sampScale + rpc_model.sampOffset;
    
end

function id = auxCompute(coeffs, L, P, H)
    id = coeffs(1)       + coeffs( 6).*L.*H + coeffs(11).*P.*L.*H + coeffs(16).*P.*P.*P ...
        + coeffs(2).*L    + coeffs( 7).*P.*H + coeffs(12).*L.*L.*L + coeffs(17).*P.*H.*H ...
        + coeffs(3).*P    + coeffs( 8).*L.*L + coeffs(13).*L.*P.*P + coeffs(18).*L.*L.*H ...
        + coeffs(4).*H    + coeffs( 9).*P.*P + coeffs(14).*L.*H.*H + coeffs(19).*P.*P.*H ...
        + coeffs(5).*L.*P + coeffs(10).*H.*H + coeffs(15).*L.*L.*P + coeffs(20).*H.*H.*H;
end

    